#include <iostream>
#include <fstream>
#include <vector>

#ifdef TIEMPOS
#include <sys/time.h>
#endif

#define forn(i,n) 	 for(int i=0; i<n; i++)
#define fornd(i,n) 	 for(int i=n; i>0; i--)
#define fornx(i,x,n) for(int i=x; i<n; i++)

using namespace std;

// Si la matriz que recibe strassen como parametro es de n x n,
// si n < CASO_BASE_STRASSEN entonces se multiplica con el algoritmo
// trivial. Default = 129
int CASO_BASE_STRASSEN = 129;

void strassen(int **m, int **a, int **b, int start_a_x, int start_a_y, int start_b_x, int start_b_y, int n, int start_m_x, int start_m_y);
void invertir(int ** a, int n);

int main()
{

    return 0;
}

void invertir(int ** a, int n)
{
    forn(f, n)
    {
        forn(c, n)
			a[f][c] = a[f][c]==1 ? 0 : 1;
    }
}

void sumar(int **m, int **a, int **b, int start_a_x, int start_a_y, int start_b_x, int start_b_y, int n, int start_m_x, int start_m_y)
{
    forn(f, n)
    {
        forn(c, n)        
            m[start_m_x + f][start_m_y + c] = a[start_a_x+f][start_a_y+c] + b[start_b_x+f][start_b_y+c];        
    }
}


void restar(int **m, int **a, int **b, int start_a_x, int start_a_y, int start_b_x, int start_b_y, int n, int start_m_x, int start_m_y)
{
    forn(f, n)
    {
        forn(c, n)        
            m[start_m_x + f][start_m_y + c] = a[start_a_x+f][start_a_y+c] - b[start_b_x+f][start_b_y+c];        
    }
}

// Estoy suponiendo n potencia de 2
// Rango [start, start+n)
void strassen(int **m, int **a, int **b, int start_a_x, int start_a_y, int start_b_x, int start_b_y, int n, int start_m_x, int start_m_y)
{
    if (n<CASO_BASE_STRASSEN)
    {
        forn(i, n)
        {
            forn(j, n)
            {
                m[start_m_x+i][start_m_y+j] = 0;

                forn(k, n)
					m[start_m_x+i][start_m_y+j] += a[start_a_x+i][start_a_y+k] * b[start_b_x+k][start_b_y+j];
            }
        }
    }
    else
    {
        int **m1 	  	= new int * [n/2];
        int **temp_m1 	= new int * [n/2];
        int **m2 		= new int * [n/2];
        int **m3 		= new int * [n/2];
        int **m4 		= new int * [n/2];
        int **m5 		= new int * [n/2];
        int **m6 		= new int * [n/2];
        int **temp_m6 	= new int * [n/2];
        int **m7 		= new int * [n/2];
        int **temp_m7 	= new int * [n/2];

        forn(i, n/2)
        {
            m1[i] 		= new int[n/2];
            temp_m1[i] 	= new int  [n/2];
            m2[i] 		= new int[n/2];
            m3[i] 		= new int[n/2];
            m4[i] 		= new int[n/2];
            m5[i] 		= new int[n/2];
            m6[i] 		= new int[n/2];
            temp_m6[i] 	= new int[n/2];
            m7[i] 		= new int[n/2];
            temp_m7[i] 	= new int[n/2];
        }

        sumar(m1,a,a,start_a_x,start_a_y,start_a_x+n/2,start_a_y+n/2,n/2,0,0);
        sumar(temp_m1,b,b,start_b_x,start_b_y,start_b_x+n/2,start_b_y+n/2,n/2,0,0);
        strassen(m1,m1,temp_m1,0,0,0,0,n/2,0,0);

        sumar(m2,a,a,start_a_x+n/2,start_a_y,start_a_x+n/2,start_a_y+n/2,n/2,0,0);
        strassen(m2,m2,b,0,0,start_b_x,start_b_y,n/2,0,0);

        restar(m3,b,b,start_b_x,start_b_y+n/2,start_b_x+n/2,start_b_y+n/2,n/2,0,0);
        strassen(m3,a,m3,start_a_x,start_a_y,0,0,n/2,0,0);

        restar(m4,b,b,start_b_x+n/2,start_b_y,start_b_x,start_b_y,n/2,0,0);
        strassen(m4,a,m4,start_a_x+n/2,start_a_y+n/2,0,0,n/2,0,0);

        sumar(m5,a,a,start_a_x,start_a_y,start_a_x,start_a_y+n/2,n/2,0,0);
        strassen(m5,m5,b,0,0,start_b_x+n/2,start_b_y+n/2,n/2,0,0);

        restar(m6,a,a,start_a_x+n/2,start_a_y,start_a_x,start_a_y,n/2,0,0);
        sumar(temp_m6,b,b,start_b_x,start_b_y,start_b_x,start_b_y+n/2,n/2,0,0);
        strassen(m6,m6,temp_m6,0,0,0,0,n/2,0,0);

        restar(m7,a,a,start_a_x,start_a_y+n/2,start_a_x+n/2,start_a_y+n/2,n/2,0,0);
        sumar(temp_m7,b,b,start_b_x+n/2,start_b_y,start_b_x+n/2,start_b_y+n/2,n/2,0,0);
        strassen(m7,m7,temp_m7,0,0,0,0,n/2,0,0);

        sumar(m,m1,m4,0,0,0,0,n/2,0,0);
        restar(m,m,m5,0,0,0,0,n/2,0,0);
        sumar(m,m,m7,0,0,0,0,n/2,0,0);

        sumar(m,m3,m5,0,0,0,0,n/2,0,n/2);

        sumar(m,m2,m4,0,0,0,0,n/2,n/2,0);

        restar(m,m1,m2,0,0,0,0,n/2,n/2,n/2);
        sumar(m,m,m3,n/2,n/2,0,0,n/2,n/2,n/2);
        sumar(m,m,m6,n/2,n/2,0,0,n/2,n/2,n/2);

        forn(i, n/2)
        {
            delete [] m1[i];
            delete [] temp_m1[i];
            delete [] m2[i];
            delete [] m3[i];
            delete [] m4[i];
            delete [] m5[i];
            delete [] m6[i];
            delete [] temp_m6[i];
            delete [] m7[i];
            delete [] temp_m7[i];
        }

        delete [] m1;
        delete [] temp_m1;
        delete [] m2;
        delete [] m3;
        delete [] m4;
        delete [] m5;
        delete [] m6;
        delete [] temp_m6;
        delete [] m7;
        delete [] temp_m7;

    }
}

void mul(vvdouble& A, vvdouble& B, vvdouble& res)
{
	forn(i, n)
	{
		forn(j, n)
		{
			double bla = 0;
			forn(k, n)
			{
				bla += A[i][k] * B[k][j];
			}
			res[i][j] = bla;
		}
	}
}

void intercambiarFilas(vvdouble& A, vdouble&  v, vvdouble& res )
{
	forn(i, n)
	{
		res[i] = A[v[i]];
	}
}
void invertirPLU(vvdouble& A, vvdouble& invA)
{
	vdouble P;
	vvdouble LU;
	vvdouble invU;
	vvdouble invL;
	vvdouble auxA;

	descomposicionPLU(A, P, LU);
	invertirTriangSup(LU , invU);
	invertirTriangInf(LU , invL);
	mul (invU, invL, auxA);
	intercambiarFilas(auxA, P, invA);
	return;
}


/*Primer tiro: "Disparamos" una matriz inversible con un K feo ...
 * disparo es Ax 
 * Tiros siguietes:
 * se levanta la informacion de la posicion del otro y se analiza de tal forma que se realize un mejor disparo
 * */
